#########################################################################
# Comparative Cryopreservation in Amphibians
#########################################################################
library(multcomp)
stderr <- function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x)))

#Import Data
setwd("C:/Dropbox/Memphis/Comparative Cryo")
Comp_Motility<-read.csv("Comparative_CPA_Motility_2.csv", header=TRUE)
summary(Comp_Motility)

# Reordering variables
Comp_Motility$Species <- factor(Comp_Motility$Species, levels = c("ANFO","ANBA","LIPI","LISE"))

# Subsetting to create a dataframe with FRESH SPERM data only
DT510 <- subset(Comp_Motility, Comp_Motility$Cryo_treatment %in% "DT0510")

#################################################
# Summary Stats
#################################################

summary(DT510$Species)
aggregate(DT510$Weight, by=list(DT510$Species), mean)
aggregate(DT510$Weight, by=list(DT510$Species), stderr)

aggregate(DT510$Total_fresh_motility, by=list(DT510$Family), mean)
aggregate(DT510$Fresh_moving_forward, by=list(DT510$Family), mean)
aggregate(DT510$Fresh_Concentration, by=list(DT510$Species), mean)

summary(Comp_Motility)
aggregate(Comp_Motility$PercentMotRemaining, by=list(Comp_Motility$Species), mean)
aggregate(Comp_Motility$PercentForwardRemaining, by=list(Comp_Motility$Species), mean)


##############################################################################
# Data analyses
##############################################################################


# Test for normality
shapiro.test(DT510$Total_fresh_motility)
shapiro.test(DT510$Fresh_moving_forward)
shapiro.test(DT510$Fresh_Concentration)
shapiro.test(Comp_Motility$PercentMotRemaining)
shapiro.test(Comp_Motility$PercentForwardRemaining)

# Test for homogeneity of variances
bartlett.test(DT510$Total_fresh_motility~DT510$Status)
bartlett.test(DT510$Total_fresh_motility~DT510$Family)

bartlett.test(DT510$Fresh_moving_forward~DT510$Status)
bartlett.test(DT510$Fresh_moving_forward~DT510$Family)

bartlett.test(DT510$Fresh_Concentration~DT510$Status)
bartlett.test(DT510$Fresh_Concentration~DT510$Family)

bartlett.test(Comp_Motility$PercentMotRemaining~Comp_Motility$Status)
bartlett.test(Comp_Motility$PercentMotRemaining~Comp_Motility$Family)
bartlett.test(Comp_Motility$PercentMotRemaining~Comp_Motility$Cryo_treatment)

bartlett.test(Comp_Motility$PercentForwardRemaining~Comp_Motility$Status)
bartlett.test(Comp_Motility$PercentForwardRemaining~Comp_Motility$Family)
bartlett.test(Comp_Motility$PercentMotRemaining~Comp_Motility$Cryo_treatment)


##############################################################
# FRESH 
##############################################################
################################
# Fresh Motility ~ Status + Family
################################
# Using GLM quasibinomial for proportional data and corrected for over-dispersion
# Using DT5-10 subset, so that there is only 1 value per individual

df2<-glm(cbind(Total_fresh_motility,100-Total_fresh_motility)~Status+Family,
         data=DT510,
         family=quasibinomial)
summary(df2)
anova(df2, test = "F")
summary(glht(df2, mcp(Species="Tukey")))

# Analysis of Deviance Table
# Model: quasibinomial, link: logit
# Response: cbind(Total_fresh_motility, 100 - Total_fresh_motility)
# Terms added sequentially (first to last)
#  Df Deviance Resid. Df Resid. Dev      F    Pr(>F)    
#  NULL                       63    1016.46                     
#  Status  1  237.316        62     779.14 20.108 3.292e-05 ***
#  Family  1   37.943        61     741.20  3.215   0.07792 .  
#  ---
#  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



################################
# Fresh forward motility
################################
# Using GLM quasibinomial for proportional data and corrected for over-dispersion
# Using DT5-10 subset, so that there is only 1 value per individual

ff1<-glm(cbind(Fresh_moving_forward, 100-Fresh_moving_forward)~Status+Family,
         data=DT510,
         family=quasibinomial)
summary(ff1)
anova(ff1, test = "F")
summary(glht(ff1, mcp(Species="Tukey")))


# Analysis of Deviance Table
# Model: quasibinomial, link: logit
# Response: cbind(Fresh_moving_forward, 100 - Fresh_moving_forward)
# Terms added sequentially (first to last)
# Df Deviance Resid. Df Resid. Dev      F    Pr(>F)    
# NULL                       63     2471.6                     
# Status  1  1009.56        62     1462.1 53.233 7.172e-10 ***
# Family  1   191.05        61     1271.0 10.074  0.002357 ** 
#   ---
#   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


################################
# Concentration
################################
# Using GLM quasiPoission for count data and corrected for over-dispersion
# Using DT5-10 subset, so that there is only 1 value per individual

con1<-glm(Fresh_Concentration~Status+Family,
          data=DT510,
          family=quasipoisson)
summary(con1)
anova(con1, test = "F")
summary(glht(con1, mcp(Species="Tukey")))


# Analysis of Deviance Table
# Model: quasipoisson, link: log
# Response: Fresh_Concentration
# Terms added sequentially (first to last)
# Df  Deviance Resid. Df Resid. Dev      F    Pr(>F)    
# NULL                        63  668547266                     
# Status  1  85820372        62  582726894 10.017   0.00242 ** 
# Family  1 160221455        61  422505439 18.701 5.755e-05 ***
#   ---
#   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1






##############################################################
# Recovery Rate 
##############################################################
################################
# Recovery Rate for Motility 
################################
ChangeMot1<-glm(cbind(PercentMotRemaining, 100-PercentMotRemaining)~
                  Status+Family+Cryo_treatment,
                data=Comp_Motility,
                family=quasibinomial)
summary(ChangeMot1)
anova(ChangeMot1, test = "F")
summary(glht(ChangeMot1, mcp(Cryo_treatment="Tukey")))

# Analysis of Deviance Table
# Model: quasibinomial, link: logit
# Response: cbind(PercentMotRemaining, 100 - PercentMotRemaining)
# Terms added sequentially (first to last)
# Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
# NULL                             191     7332.6                      
# Status          1     28.0       190     7298.5  1.1778 0.2792    
# Family          1    121.4       189     7177.1  5.1055 0.0250 *  
# Cryo_treatment  2   4014.1       187     3163.0 84.4159 <2e-16 ***
#  ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Simultaneous Tests for General Linear Hypotheses
# Multiple Comparisons of Means: Tukey Contrasts
# Fit: glm(formula = cbind(PercentMotRemaining, 100 - PercentMotRemaining) ~ 
#            Status + Family + Cryo_treatment, family = quasibinomial, 
#           data = Comp_Motility)
# Linear Hypotheses:
#   Estimate Std. Error z value Pr(>|z|)    
#   DT1010 - DT0510 == 0  -2.4754     0.3119  -7.937   <0.001 ***
#   DT1510 - DT0510 == 0  -3.7301     0.5272  -7.076   <0.001 ***
#   DT1510 - DT1010 == 0  -1.2548     0.5843  -2.148   0.0752 .  
# ---
#   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#   (Adjusted p values reported -- single-step method)



################################
# Recovery Rate for Forward Progressive Motility 
################################
fwd1<-glm(cbind(PercentForwardRemaining, 100-PercentForwardRemaining)~
            Status+Family+Cryo_treatment,
          data=Comp_Motility, 
          family=quasibinomial)
summary(fwd1)
anova(fwd1, test="F")
summary(glht(fwd1, mcp(Cryo_treatment="Tukey")))

# Analysis of Deviance Table
# Model: quasibinomial, link: logit
# Response: cbind(PercentForwardRemaining, 100 - PercentForwardRemaining)
# Terms added sequentially (first to last)
# Df Deviance Resid. Df Resid. Dev       F   Pr(>F)    
# NULL                             179     7998.0                     
# Status          1     63.1       178     7934.8  0.9907   0.32095    
# Family          1    388.0       177     7546.8  6.0876   0.01458 *  
# Cryo_treatment  2   3437.6       175     4109.2 26.9672 6.182e-11 *** 
#   ---
#   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


# Simultaneous Tests for General Linear Hypotheses
# Multiple Comparisons of Means: Tukey Contrasts
# Fit: glm(formula = cbind(PercentForwardRemaining, 100 - PercentForwardRemaining) ~ 
#            Status + Family + Cryo_treatment, family = quasibinomial, 
#          data = Comp_Motility)
# Linear Hypotheses:
#   Estimate Std. Error z value Pr(>|z|)    
#   DT1010 - DT0510 == 0  -2.9942     0.6925  -4.324   <1e-04 ***
#   DT1510 - DT0510 == 0  -3.4346     0.8344  -4.117   <1e-04 ***
#   DT1510 - DT1010 == 0  -0.4405     1.0281  -0.428    0.901    
# ---
#   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#  (Adjusted p values reported -- single-step method)


